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Abstract 

A recently introduced particle-based model for fluid dynamics with continuous ve- 
locities is generalized to model fluids with excluded volume effects. This is achieved 
through the use of biased stochastic multi-particle collisions which depend on local 
velocities and densities and conserve momentum and kinetic energy. The equation 
of state is derived and criteria for the correct choice of collision probabilities are 
discussed. In particular, it is shown how a naive implementation can lead to incon- 
sistent density fluctuations. 
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1 Introduction 

Simulation studies of the structure and dynamic properties of complex liq- 
uids are often complicated by the fact that typical energy scales are on the 
order of the thermal energy and the characteristic structural length scales 
are in the range of nanometers to micrometers. The resulting large number 
of degrees of freedom and disparate length and time scales require the use of 
"mesoscale" simulation techniques which achieve high computational efficiency 
by "averaging out" irrelevant microscopic details while retaining the essential 
features of the microscopic physics on the length scales of interest. The fact 
that the properties of these systems are strongly influenced by a delicate in- 
terplay between thermal fluctuations, hydrodynamic interactions, and possible 



Preprint submitted to Elsevier Science 



2 February 2008 



w 



Fig. 1. Collision rules. Four distinct collisions are considered: a) horizontally along 
tri, b) vertically along 0-2, c) diagonally and d) off-diagonally along 0-3 and <x 4 . w 
and Wd denote the probabilities of choosing collisions a), b) and c), d) respectively. 

Fig. 2. P n times r as a function of M 1 / 2 , measured using the microscopic stress 
tensor. Symbols show data in the range r = 0.05, . . . ,2.00. Parameters: L/a = 32, 
&bT = 1.0. The inset shows the nonideal contribution to the pressure as a function 
of (ksT) 1 / 2 . Parameters: L/a = 32, M = 10, r = 0.40. In both figures, the solid 
lines are the theoretical prediction given by Eq. (1). 

spatio-temporally varying forces, places additional stringent requirements on 
the simulation protocol. A recently introduced particle-based simulation tech- 
nique [1] — often called stochastic rotation dynamics (SRD) [2,3,4,5,6], multi- 
particle collision dynamics [7], or real-coded lattice gas [8] — is a promising 
algorithm for mesoscale simulations of this type. SRD solves the hydrody- 
namic equations of motion by following the path of fluid point-particles in 
discrete time and continuous space. Efficient multi-particle collisions which 
explicitly conserve momentum and energy enable simulations in the micro- 
canonical ensemble, while fully incorporating both thermal fluctuations and 
hydrodynamic interactions. Furthermore, its simplicity has made it possible 
to obtain accurate analytic expressions for the transport coefficients which are 
valid for both large and small mean free paths. 

The original SRD algorithm — in which collisions consist of a stochastic ro- 
tation of the relative velocities of particles in the collision cells — describes 
a fluid with an ideal gas equation of state. The fluid is therefore very com- 
pressible, and the speed of sound, c s , is low. However, this collision rule is 
not unique and other choices of collision rules can lead to nonideal behav- 
ior. In this paper we consider one such collision rule and discuss in detail 
the conditions which must be fulfilled in order to guarantee thermodynamic 
consistency. There are, however, several subtleties associated with collision 
rules of this type, such as phase space contraction and semi-detailed balance 
in this reduced phase space, which are not addressed here. In systems with 
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explicit interparticle potentials — such as the hard sphere fluid — this behavior 
can be analyzed in considerable detail. However, to the best of our knowledge, 
there are no non-trivial models with multi-particle interactions of the type 
we consider for which this has been done. The thermodynamic consistency of 
such models is not ensured a priori, and the aim of this paper is to provide 
guidelines for the construction of consistent models and to point out possible 
pitfalls. 

A more realistic modeling of dense gases and liquids can be achieved by in- 
troducing generalized excluded volume interactions between the fluid parti- 
cles. The resulting algorithm can be thought of as a coarse-grained multi- 
particle collision generalization of a hard sphere fluid, since, just as for hard 
spheres, the kinetic energy is conserved. There is no potential energy, so 
that the internal energy is the same as that of an ideal gas. Thermody- 
namic consistency therefore requires that c v = T ds/dT\ p = dk B /2, where 
d is the spatial dimension. It follows that the nonideal contribution to the 
entropy density, s n , can only depend on the density, p, so that the free en- 
ergy density f(T,p) = fideai(T,p) + Ts n (p). The equation of state is therefore 
P = pdf/dp\ T - f = Pideai + T[pds n (p)/dp - s„], and P - P ideai is strictly 
proportional to the temperature T. 



2 Model 



As in the original SRD algorithm, the solvent is modeled by N of point- 
like particles of mass m which move in continuous space with a continuous 
distribution of velocities. The system is coarse-grained into (L/a) d cells of 
a d-dimensional cubic lattice of linear dimension L and lattice constant a. 
The algorithm consists of individual streaming and collision steps. In the free- 
steaming step, the coordinates, r^t), of the solvent particles at time t are 
updated according to r«(t + r) = rj(t) + rvj(t), where Vj(t) is the velocity of 
particle i at time t and r is the value of the discretized time step. In order to 
define the collision, we introduce a second grid with sides of length 2a which 
(in d = 2) groups four adjacent cells into one "supercell". For simplicity, 
we restricted ourselves to two dimensions, but the algorithm can be easily 
extended to three dimensions. 

As proposed in Ref. [2], a random shift of the particle coordinates before 
the collision step is required to ensure Galilean invariance. All particles are 
therefore shifted by the same random vector with components in the interval 
[—a, a] before the collision step. Particles are then shifted back by the same 
amount after the collision. To initiate a collision, pairs of cells in every supercell 
are randomly selected. As shown in Fig. 1, three distinct choices are possible: 
a) horizontal (<Ti), b) vertical (cr 2 ), and c) diagonal collisions (cr 3 and cr 4 ). 
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In every cell, we define the mean particle velocity, u n = (1/M n ) J2i=i v «, 
where the sum runs over all particles, M n , in the cell with index n. The 
projection of the difference of the mean velocities of the selected cell-pairs 
on (Tj, Au = <jj • (ui — u 2 ), is then used to determine the probability of a 
collision. If Au < 0, no collision will be performed. For positive Au, a collision 
will occur with an acceptance probability which depends, in principle, on Au 
and the number of particles in the two cells, M 1 and M 2 . Indeed, the choice of 
acceptance probability, pa, determines both the equation of state and values of 
the transport coefficients, and the requirement of thermodynamic consistency 
imposes severe restrictions on the choice of pa- In Ref. [9] it was shown that 
the choice p A (M u M 2 , Am) = Q(Au) tanh(A), with A = AAuM l M 2 , where 6 
is the unit step function and A is a (small) constant leads to an equation of 
state of the required form. In this paper we explore the consequences of the 
simpler choice pa = O(A-u), which is identical to the limit A — > oo of the 
model presented in Ref. [9] . 

The collisions should conserve the total momentum and kinetic energy of the 
cell-pairs participating in the collision, and in analogy to the hard-sphere liq- 
uid, they should primarily transfer the component of the momentum which is 
parallel to the connecting vector crj. The rule we have chosen is to exchange the 
parallel component of the mean velocities of the two cells, which is equivalent 
to a "reflection" of the relative velocities [9]. The perpendicular component 
remains unchanged. Because of x — y symmetry, the probabilities for choosing 
cell pairs in the x— and y— directions are equal, and will be denoted by w. 
The probability for choosing diagonal pairs is given by Wd = 1 — 2w. w and 
Wd must be chosen to that the hydrodynamic equations are isotropic and do 
not depend on the orientation of the underlying grid. As shown in Ref. [9], 
this can be achieved only if Wd — 1/2 and w — 1/4. 



3 Equation of state and the structure factor 

The pressure can be calculated using the method described in [9]. For pa = 
0(Aw), one finds 

P = Pideal + Pn= pk B T + (6/t) {pk^T (1) 

in the limit of large M, where p = M/a 2 is the particle density and b = 
(1/4 + l/2v / 2)/2 A /7T. The first term in Eq. (1) is the ideal gas contribution 
and the second is the contribution from the collisions, the nonideal pressure, 
P n . It can be shown that P n ~ p 2 in the limit of small M. Simulation results 
for P n obtained by averaging the diagonal part of the microscopic stress tensor 
[9] were found to be in good agreement with Eq. (1), see Fig. 2. This equation 
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Fig. 3. Speed of sound as a function of r. (•) and (O) are results for k = (1,0) and 
(0, 1), respectively. The solid line is the theoretical prediction given by Eq. (2). The 
inset shows the dynamical structure factor as a function of uj for k = (2, 0), r = 0.2. 
Parameters: L/a = 128, M = 5, k B T = 1.0. 

Fig. 4. = 0)/p as a function of r (•). (O) are results obtained by numerically 

evaluating the derivative of the pressure measured using the microscopic stress ten- 
sor. The solid line is a plot of Eq. (3). k is the lowest wave vector. 



of state is not consistent with the fact that the kinetic energy is conserved, 
since P — Pideai ~ VT instead of T. Using Eq. (1), the adiabatic speed of 
sound in d = 2 is 
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We have performed simulations to determine the dynamical structure factor, 
S(k,u) = {p{k ) ui)p{k,u)) . Measuring the structure factor can be tricky and 
details are discussed elsewhere. S(k ) uo) is plotted as a function of uj in the 
inset to Fig. 3. The position of the finite frequency peaks in S(k, uj) gives the 
speed of sound. The solid vertical lines in the figure show the theoretically 
predicted positions of the frequencies for the speed of sound given in Eq. (2). 
The dashed lines show the peak positions for an ideal gas. Results for the speed 
of sound for various /c-values is shown in Fig. 3. The agreement with Eq. (2) 
is satisfactory. We have also checked that the sound speed is isotropic for the 
model. Thermodynamics provides a relation between density fluctuations and 
the derivative of the pressure; i.e. using Eq. (1), 



S(k,t = 0) = pk B T^ 5 



1 + 6/(2rVp^T) 



(3) 
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Fig. 4 compares simulation data for S(k,t = 0)/p, results for the second ex- 
pression in Eq. (3) obtained by taking the numerical derivative of the measured 
pressure, and the analytical result, the last term in Eq. (3). Data obtained by 
measuring the fluctuations in the density, S(k,0), are clearly not consistent 
with the results based on measurements of the expectation value of the di- 
agonal part of the microscopic stress tensor. As already discussed, the reason 
for this is that the equation of state is not consistent with the fact that the 
algorithm conserves kinetic energy. Note that for = Q(Au), the collision 
probability does not depend on the density, which is clearly unphysical. More- 
over, the presence of a ^fp term in Eq.(l) is strange since the first term in 
a typical virial expansion would be proportional to p 2 . The consequences of 
choosing collision rules which violate thermodynamic consistency are there- 
fore quite dramatic and easy to detect. However, as shown in Ref. [9], for the 
correct choice of collision probabilities, the algorithm can be made thermody- 
namically consistent. 
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